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In this article, we numerically solve an equation modeling the evolution of the density 
of glioma in the brain—the most malignant form of brain tumor quantified in terms of 
net rates of proliferation and invasion. We employ a non-linear heterogeneous diffusion 
logistic density model. This model assumes that glioma cell invasion throughout the brain 
is a reaction-diffusion process and that the coefficient of diffusion can vary according 
to the gray and white matter composition of the brain at that location. The analysis 
provided in this article demonstrates that using the correct finite difference scheme can 
overcome the stability issues caused by the discontinuities of the diffusion coefficient. 
We also observe that at the steady-state these discontinuities vanish. To visualize and 
investigate numerically the behavior of the evolution of tumor concentration of the glioma, 
we calculated and plotted the number of tumor cells, the average mean radial distance, and 
the speed of the tumor cells along with charting the effects of net dispersal rate and net 
proliferation rate terms versus time for different center position values of Gaussian initial 
profile for each zone (gray and white matter tissues). We have proposed two numerical 
methods, the implicit backward Euler and the averaging in time and forward differences in 
space (the Crank-Nicolson scheme), both in combination with Newton's method for solving 
the governing equations. These methods are compared in terms of their performance in 
varying time-step and mesh-discretization. The Crank-Nicolson implicit method is shown 
to be the better choice to solve the equation. 

© 2015 Elsevier Ltd. All rights reserved. 


1. Introduction 

Cancer is a complex disease which leads to the uncontrolled growth of abnormal cells, destruction of normal tissues 
and invasion of vital organs. Uncontrolled evolution of cells leads to tumor formation but the tumors can be benign and 
malignant. The distinction between benign and malignant tumors is based on many criteria: degree of differentiation, rate 
of growth, local invasion and metastatic ability (Prayson [1]). 

Mathematical modeling of tumor growth by diffusion has been studied by many researchers (Clatz et al. [2], Cruywagen 
et al. [3], Glass [4], Rockne et al. [5], Tracqui et al. [6], and Woodward et al. [7]). Unfortunately these works did not consider 
the spatial heterogeneity of the brain tissue in terms of gray and white matter. Giese et al. [8] experimentally established 
the result that tumor cells moved faster in white matter than in gray matter. Swanson et al. [9] took into account the spatial 
heterogeneity of the brain tissue by defining the diffusion coefficient D as a function of the spatial variable differentiating 
regions of gray and white matter and observed that migration in white brain tissue was faster than in gray tissue. 
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In a reaction-diffusion model, diffusion and proliferation are the main biological behaviors. This approach models tumor 
cell infiltration into the surrounding brain tissue. The proliferation is a function representing a reactive behavior and it sim¬ 
ulates tumor cell growth and death (Murray [10]). In the literature, the problem of brain tumor growth has been formulated 
as a reaction-diffusion process, in which the rate of change of tumor cells is given by the net proliferation and diffusion 
of tumor cells (Bacaer [11], Drasdo [12], Hatzikirou et al. [13], Harpold [14], Mansuri [15], Murray [16], Swanson [17,18]). 
Swanson’s thesis [17] is especially important because it encompasses both the mathematical and medical issues at work. 
Reaction-diffusion equations with logistic nonlinearity were introduced in the pioneering works of Fisher [19], and Kol¬ 
mogorov, Petrovskii and Piskunov [20], The latter model is known as the KPP equation. 

Different kinds of cancerous growth such as exponential, logistic and Gompertz are detailed in Murray [16]. Population 
growth equations are commonly used for the proliferation rate/(c) as summarized either as the exponential law 

f(c) = pc(t,x) 


or the Verhulst (or logistic) law 

c(t,x) 


c(t,x) 


/(c) = pc(t,x) |^1 - 
or Gompertz law 

/(c) = -pc(t, x) In 


e k/d 


where p is the net proliferation rate in units day -1 , c max represents the carrying capacity of the tissue, providing an upper 
limit on the number of tumor cells capable of occupying any cubic millimeter of brain, k is the growth rate of tumor and d 
is a density coefficient, and e k ^ d is the carrying capacity. Here, k gives an exponential increase when the tumor size c(t , x) is 
small and d, the decay constant, damps out the growth rate when c(t, x) is large. 

The exponential growth is the simplest proliferation law and is suited for quantifying the growth of small tumors during 
a short period. Gompertz [21] showed that the growth, which was initially exponential, was later limited to an asymptotic 
rate. In the Gompertzian mode, tumors develop with a growth rate that decreases with the tumoral growth. 

The estimation of parameters such as cell proliferation (p), diffusion coefficient D, and carrying capacity c max is important 
in understanding the behavior of glioma. Since medical imaging techniques can detect regions containing tumor cells, but 
only if the number of cells are above some threshold (Konukoglu [22]), net rates of proliferation (p) and dispersal (diffusion 
coefficient) D have been estimated from features of pretreatment magnetic resonance (MR) images to predict tumor growth 
(Tracqui et al. [6], Swanson et al. [9], Powathil et al. [23], Swanson [18], Szeto et al. [24], Konukoglu et al. [25], Roniotis 
et al. [26]). 

The two general numerical methods used to solve the reaction-diffusion equation in the context of tumor growth are 
finite elements and finite differences. Clatz et al. [27], Hogea et al. [28], and Rockne et al. [29] used finite-element methods 
to solve anisotropic diffusion equation. But, the finite difference scheme is easier to implement and the pixel structure of 
digital images provides a natural regular grid. However, the anisotropic diffusion equation is a complex second order partial 
differential equation (PDE) and simple finite difference is not appropriate. This second order PDE was reduced to simpler 
first order PDEs by using the concept of manifolds (Cobzas et al. [30] and Konukoglu et al. [25]). But the solution of the 
equation gives us tumor delineation area instead of the tumor cell density. Jbabdi et al. [31] used chain rule expansion 
to solve anisotropic diffusion equation, however, Mosayebi et al. [32] proved that chain rule expansion was unstable. In 
Roniotis [33], anisotropic heterogeneous diffusion case study, the abstract model used was non-linear. However, for their 
computation, the net cell proliferation rate function was zero which effectively rendered the model linear. They made a 
detailed error analysis in this linear case. Lolas [34] used finite difference methods to solve isotropic diffusion equation. An 
isotropic diffusion equation was first proposed by Swanson [17] in her Ph.D. thesis assuming the net growth of glioma cells 
was linear. Later, Swanson et al. [35], Roniotis et al. [26], Papadomanolaki and Saridakis [36] utilized the finite difference 
methods to solve heterogeneous version of this diffusion equation. 

We should note that tumor environment is combination of many constituents, such as, different types of tumor cells, 
blood vessels, nutrients, immune system cells, healthy cells, among others. There are numerous researchers modeling the 
tumor evolution using complex models that have 4-, 6-, and 10-species. For instance, the model developed by Cristini 
et al. [37] depicted the evolution of the nutrients as a reaction-diffusion equation and delineated the mixture as being 
composed of tumor cells and water. As the model did not consider the dead cells as a separate constituent, when the cells 
inside the tumor died because of lack of oxygen, water occupied their original position, and the characteristic necrotic core 
was not represented. The model developed by Wise et al. [38] splits the tumor into viable and dead cells. This modification 
allowed the representation of cells dying inside the tumor due to lack of nutrients and oxygen. 

The model developed by Hawkins-Daarud et al. [39] was a four-species tumor growth model, in which the extracellular 
water was divided into nutrient rich and nutrient poor parts. However, it did not reproduce the necrotic core inside the 
tumor. A universal hybrid 10-species tumor growth model of this general type was recently developed by Lima et al. [40], 
in which tumor cells were divided into proliferative, hypoxic (quiescent) and necrotic states: healthy cells, nutrient, tumor 
growth factor and endothelial cells all incorporated in the mixture. Following the work of Lima et al. [40], Lima et al. [41] 
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developed a system of stochastic partial differential equations governing a six-species tumor growth model accounted for 
random perturbations that arose from modeling uncertainties in the parameters. They considered an avascular model with 
six constituents derived directly from the hybrid 10-species model developed in Lima et al. [40] and simulated the stochastic 
behavior of cellular and macrocellular events affecting the evolution of avascular cancerous tissue using a mixed finite 
element method and a stochastic collocation scheme. 

An extensive review of diffusive modeling of glioma evolution can be found in Roniotis et al. [26], Hatzikirou et al. [13], 
Giatili and Stamatakos [42], Belmonte-Beitia et al. [43], Branco et al. [44] and references therein. 

Recently, in 2014, a book titled by ‘Computational Surgery and Dual Training: Computing, Robotics and Imaging’ and 
edited by Garbey et al. [45] was published, and in the Chapter 20 of this book, Sodt et al. [46] employed a finite volume 
and backward Euler method to discretize the reaction-diffusion equation in 3D. They quantitatively compared simulated 
and actual 3D tumor volumes to assess the impact of anisotropic diffusion on the accuracy of their model. By means of a 
mathematical model for glioma growth, they saw that anisotropic migration of tumor cells along white matter fiber tracts 
could influence the overall diffusive growth pattern of gliomas. 

In this study, our aim is not to discredit the multi species tumor growth models mentioned above but to focus on the 
methodology of handling stability issues due to the discontinuous diffusion coefficient D(x) in a non-linear heterogeneous 
diffusion logistic density model. 

Our model has been extensively studied (Swanson [17], Swanson et al. [9], and Swanson [47]) and found to be very 
effective in simulating the behavior of real malignant brain tumors in the time frame for gliomas. The model we are 
considering in the present work is an analog of this 3-dimensional problem (Sodt et al. [46]) in one spatial dimension. 

In this paper, the effects of radiotherapy and chemotherapy are neglected. We adopt a mathematical density model 
of glioma which assumes that glioma cell invasion throughout the brain is a diffusion process and that the coefficient of 
diffusion can vary in space depending on the gray and white composition of the brain at that location. In previous studies 
(Swanson [17,35]) linear growth was considered. Later on, in 2008, Swanson et al. [48] considered logistic growth but with 
constant diffusion coefficient. We use a non-linear heterogeneous diffusion logistic density model. Unlike Swanson [48], our 
model is 1-dimensional, however, we assume that the diffusion coefficient D depends on the tissue environment to reflect 
the observation that glioma cells exhibit higher motility in white matter than in gray matter, thus we consider D as a function 
of x having a finite number of discontinuities and this results a lack of smoothness in the solution. This is an obstacle for 
standard numerical procedures, see Section 3.1. However, at the steady-state, these discontinuities vanish. 

The formulation and the numerical procedure are described in Sections 2 and 3, respectively. An important question is 
how the diffusion term is discretized. This is addressed in Section 3. 1. The numerical simulations are presented in Section 4.1, 
the stability condition is derived in Section 4.2, and the analysis is done in Section 4.3. Finally, the conclusion is presented 
in Section 5. 


2. Formulation of the problem 


Using parameter ranges obtained from Tracqui et al. [ 6 ], along with the evidence that suggests glioma cells migrate with 
greater velocity through white matter of the brain than through gray (Giese et al. [8]), we extend the model used by Swanson 
et al. [48] by assuming that the diffusion coefficient D depends on the tissue environment representing the net motility of 
tumor cells or the Fickian random dispersal rate of glioma cells in undifferentiated brain tissue: 


3c(t, x) 
3 1 


3 

3x 



3c(f, x)\ 
3x ) 


+ pc{t,x) 


c(t,x) 

Qnax 


and 


D(x) = 


D g = 0.13 mm 2 /day 
D w = 0.65 mm 2 /day 
D g =0.13 mm 2 /day 


0 < x < 7.5 (gray region) 

7.5 < x < 42.5 (white region) 

42.5 < x < 50 (gray region) 


(see Fig. 1(b)) with the Neumann boundary conditions (zero-flux) 


( 2 . 1 ) 


c*(0, t) = 0 
c x (50, t) = 0 

and the initial condition is taken as Gaussian initial tumor profile (see Fig. 1(a)) 


( 2 . 2 ) 


1 i 

c(0,x)=g(x)= _ e 2 V e / , x e (0, 50) 

Vzjre 


(2.3) 


where c = c(t, x) denoted the tumor concentration of glioma cells at time t and spatial location x and x 0 = 25 mm (the 
middle of the considered interval) and e = 0.01 suggested by Becker et al. [49]. Here, p is the net proliferation rate in units 
day -1 and Fickian diffusion has been used to quantify the random motility of a variety of invading cells (Cozens-Roberts 
et al. [50]). A factor of 5, D w = 5 D g , was used by Swanson et al. [5 1 ] but this may vary from patient to patient. The (maximum) 
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Fig.l. (a) The initial distribution c(0, x) =-^=-e 2 V £ > versus x with x 0 = 25 mm (the middle of the considered interval), e = 0.01, and Ax = 0.5 mm. 
(b) The diffusion coefficient D(x) versus x with Ax = 0.5 mm. 


net proliferation rate p includes both birth and death rates and assumes logistic growth with a tissue-carrying capacity c max . 
The net proliferation term includes mitosis, apoptosis, and other cell loss mechanisms. There is a spatial heterogeneity in 
the net proliferation term resulting from the effect of c max (Szeto [24]). The model assumes logistic growth of the tumor cell 
population, so that the net proliferation rate p is lower in regions of high cell density (where c ~ c max ) than in regions of 
low cell density (where c c max ). As Swanson [17], Harpold et al. [14], Swanson and Alvord [52] suggested, this model 
served well in capturing tumor dynamics in vivo thus far in many iterative comparisons of theory and the real world. This 
model has proven to be an accurate predictor of simulated tumor growth with respect to anatomical changes that can be 
imaged by Tl-weighted gadolinium enhanced (TIGd) and T2-weighted MRI (Swanson et al. [47]). 

Eq. (2.1) states that tumor cells at a given x either originate from moving to a close location by passive diffusion or from 
cell division by means of a logistical proliferation model. Eq. (2.1) is a 1-dimensional parabolic partial differential equation. 
Theg function can be considered as the distribution of the initial cell concentration. If the function g is continuous on [0, 50], 
then the boundary value problem (2.1) has a unique solution. 

This model corresponds to the following law: the rate of change of the tumor cell population density is equal to the 
diffusion of the tumor cells in white and gray matter plus net proliferation of the tumor cells, (Swanson et al. [9]). This PDE 
(2.1) describes the macroscopic growth and invasion of gliomas under a variety of treatment conditions (Rockne et al. [29], 
Swanson et al. [53]). Microscopic growth rates describe cell divisions and invasions by means of interactions between tumor 
cells and their surrounding tissue (Cristini et al. [54], Patel et al. [55]). 

Swanson [ 17] in her Ph.D. thesis used the model given in Eq. (2.1) in order to more accurately reflect the spatial limitation 
of cellular proliferation and inherent heterogeneity in the brain by allowing the diffusion coefficient D to depend on the tissue 
environment, D = D(x), however she favored exponential growth for the proliferation rate term/(c) = p c(t, x). The case 
of exponential growth/(c) = p c(t, x) and constant diffusion D was studied (Cruywagen et al. [3], Tracqui et al. [6], Cook 
et al. [56], Woodward et al. [7], Swanson et al. [9,47], Clatz et al. [27]). 

In 2008, Swanson et al. [48] considered logistic growth but with constant diffusion coefficient. Here, we assume that the 
diffusion coefficient D depends on the tissue environment, namely gray and white matters, thus D is considered as a function 
of x having finite number of discontinuities and this results in the lack of smoothness in our model solution. This creates a 
difficulty and is hard to handle with standard numerical procedures, see Section 3.1. 

Fickian diffusion has been used to quantify the random motility of a variety of invading cells (Cozens-Roberts et al. [50]) 
and is implemented in our model to represent the spatial invasion of glioma cells at a rate D(x) which varies depending on 
the location in the brain to allow for an increased velocity of migration through white matter compared with gray matter. 
This migration rate in white matter is generally thought to be faster than that in gray matter: D w > D g (Giese et al. [8]). 

Eqs. (2.1) and (2.2) are solved by applying two different finite-difference methods, namely, implicit backward Euler 
method and the Crank-Nicolson method. 

3. Numerical procedure 


3.1. The Chain-rule discretization 


It might be tempting to apply the chain rule to rewrite 


9 

9x 



9c(t,x)\ 
9x ) 
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dD(x) dc(t, x) d 2 C{t,X ) 

dx dx 9x 2 

and then to apply the discretization (Roniotis [33]). However, this is not the best approach as mentioned in Mosayebi 
et al. [32] in which they proved for 2D and 3D cases that chain rule expansion violated the stability issues and suggested 
the directional discretization method which was both theoretically and practically stable if the condition on tensor was 
satisfied. It is better to discretize the physical problem given in Eq. (2.1) directly. This can be done by first approximating 
D(x) 3C g t x ,x) at points halfway between the grid points, using a centered approximation 

9 c(X;+i/ 2, t) ^ / Q+i — Cj 


D(x i+ 1 /2 ) 


; Dj+l/2 


dx ‘' V Ax 

and the similar approximation at Xj_i/ 2 , where Ax 

at the grid point x,- as suggested in LeVeque [57, p. 36]: 


: x i+] — X,. Subtracting these then gives a centered approximation to 




1 

Ax 

1 


D, 


1 + 1/2 


Q+i Q 
Ax 


— Di-1/2 


Cj Cj — \ 

Ax 


(Ax) : 


[D, -1/2 Q-i — (Dj- 1/2 + D, + i/2) Q + D i + i/2 C i+ i]. 


(3.4) 


3.2. The implicit backward Euler (IBE) method 

In this section, we consider an algorithm called the implicit backward Euler method, which results from substituting 
the backward difference approximation to discretize the PDE (2.1). The spatial domain [a, b] = [0, 50] (measured in 
millimeters) is divided into M sections, each of length Ax = . x* = (i — 1) A x and i = 1, ..., M — 1, and each of 

duration At = 0.01 ~ 14.4 min. or At = 0.02 » 28.8 min. The implicit backward Euler scheme is applied to Eq. (2.1) 
using Eq. (3.4), yielding 


At 



'i-l /2 cfff — (Di- 1/2 + Dj + i/ 2 ) cf +1 + D;+ 1/2 cf_+ 
2 (Ax) 2 



(3.5) 


The Neumann type of boundary conditions given in Eq. (2.2) are approximated by 


,.17+1 


- C, 


n+1 

0 


-n+1 


2 Ax 


0, 


- c; 


n+1 


Thus, we take Cq + 1 = c^ -1 " 1 and cj^ +1 = c ^_ 1 2 . Eq. (3.5) can be rewritten as 


^n+1 


2 Ax 

^n+i 


0. 


n+i 


At 


2(Ax) : 


(Di- 1/2 O — (Dj-i/2 + D, + 1/2 ) c " +1 + Dj_|_ \J2 c^ 1 ) + 


At p 


.n+1 


n n+1 


1 - 


- C, n+1 + cf 


(3.6) 


where cf = c(t n , X;) = c(n At, i A x). 

We refer this method as Method IBE. In Eq. (3.6), cf +1 , n = 0,... and i = 1, ..., M — 1 are the unknowns and cf is the 
initial data. At each time step, M — 1 non-linear algebraic equations are solved to find M — 1 unknowns for fixed values of 
D(x) and p. That system is solved by using Newton’s method. The Method IBE is first-order accurate in time and second-order 
accurate in space. 


3.3. The Crank-Nicolson (CN) scheme 

To have a second-order accurate method in both time and space, the average of the central difference approximations of 
the spatial derivatives of Eq. (2.1) at the two points n+1 and n is taken, yielding 


At 


Di- 1/2 cf_7 — (D,_i /2 + D i+ i/ 2 ) cf +l + D i+1 /2 cf+f 
2 (Ax) 2 


Di- 1/2 cf_, — (D,_i /2 + D; +1 / 2 ) cf + D i+ \/2 c" +1 
2 (Ax) 2 


.n+i 


+ cf 


1 - 


= 0 


(3.7) 
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(a) The Crank-Nicolson. 


x 

(b) The implicit backward Euler. 


Fig. 2. The values of the tumor concentration c(t, x) of glioma cells versus x: p = 0.012 day \ At = 0.01 « 14.4 min., c max = 62.5 cells/mm, 
c™ x = 39.89 cells/mm and Ax = 0.5 mm, t = 1 (asterisks), t = 5 (circles), t = 10 days (dashed), and t = 50 days (dash-dot). 


which leads to the so-called Crank-Nicolson method. We refer this method as Method CN. In Eq. (3.7), c" +1 , i = 1,..., M — 1 

are the unknowns for n = 0,_At each time step, M — 1 non-linear algebraic equations are solved to find M — 1 unknowns 

for fixed values of D(x) and p. That system is solved by using Newton’s method. The Crank-Nicolson method is second-order 
accurate in both time and space and unconditionally stable. It is recalled that an algorithm is called numerically stable if an 
error, whatever its cause, does not grow to be much larger during the calculation. 


4. Results and discussion 


4.1. Numerical simulations 


For all the numerical experiments we describe here, we used 101 spacial grids on [0, 50] with a time step At = 
0.01 day ~ 14.4 min. The CPU time for solving the partial differential equations in 1280 x 0.01 time steps is about 94 min. 
with Method CN on Linux computer which had an Intel Pentium D 3.00 GHz processor, with 3.50 GB RAM. 

Throughout simulations, a five-fold difference in the diffusion coefficients in gray and white matter is used: D u ,j,; te = 
5Dgr ay and the growth rate, p = 0.012 day -1 in the model, as suggested for high-grade gliomas by Swanson et al. [48], 

All the computational runs are made for c max = 62.5 cells/mm and for 


c(0, x) 


1 i A~*o) 


s/Itte 


where x 0 = 25 mm (the middle of the considered interval, however we have considered other values of x 0 as well, 
see Section 4.3.1) and s = 0.01. The initial distribution c(0,x) has a maximum at x = x 0 with a value of 39.89 » 
^-' o oi cells/mm, namely c™ x , i.e. we assumed that the tumor has grown to about 39.89 cells/mm as a local density before 
it begins to diffuse, see Fig. 1(a). 

However, the cases of c max = 5000 cells/mm with c™ x = 3000 cells/mm and c™ x = 500 cells/mm have not shown 
any changes in terms of the graphics but only made a difference on how many days it took to reach the carrying capacity 
and yielded 1648 days and 1721 days, respectively. 

Note that our model shows the tumor cell population reaching a maximum cell number that the tissue can carry, see 
Figs. 2, 4 and 7(a). However, in reality, after this point, new mutations will continue to occur within the nuclear DNA, 
providing advantages in their proliferation, survival and invasion. 

Fig. 2(a), (b), and Table 1 tell us that the Method CN is not only more efficient in terms of errors but also much faster in 
terms of iteration steps. We can see from Fig. 2(a) and (b) that it takes 1 day for the initial tumor cell concentration to be 
reduced to 7.24 cells/mm from 39.89 cells/mm in Method CN. However, in Method IBE this concentration is 10.62 cells/mm. 
As cells divide, their mutation pattern gives a growth advantage to some of them. This phenomena creates a new 
subpopulation of cells with higher proliferation rates. As time increases, this subpopulation of cells increases in number 
within the tumor mass, whereas those cells with lower proliferation rates decrease in number, see Figs. 2, 3 and 6 and 
Table 1. 

Fig. 3 shows the tumor concentration of glioma cells versus x using the Method CN on the left and Method IBE on the 
right for p = 0.012 day -1 , At = 0.01 ~ 14.4 min. and Ax = 0.5 mm, for the time intervals 100 < t < 200 and 
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(a) The Crank-Nicolson. (b) The implicit backward Euler. 
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c(t.x) 
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46 
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x x 

(c) The Crank-Nicolson. (d) The implicit backward Euler. 

Fig. 3. The values of the tumor concentration of glioma cells versus x: p = 0.012 day -1 , At = 0.01 « 14.4 min., c max = 62.5 cells/mm, 
c 0 max = 39.89 cells/mm and Ax = 0.5 mm. (a) and (b) t = 100 (solid), 150 (asterisks), and 200 days (plus mark), (c) and (d) t = 1050 (dotted), 
1200 (dashed), and 1280 days (dash-dot). 



Fig. 4. The evolution of the tumor concentration of glioma cells by using the Crank-Nicolson method: p = 0.012 day \ A t = 0.01 « 14.4 min., 
c max = 62.5 cells/mm, c™ ax = 39.89 cells/mm and Ax = 0.5 mm, 0 < x < 50 ,M= 100, 0 < t < 1280. 


1050 < t < 1280, respectively. From Fig. 3(a), it can be observed that by t = 200 days, the efficiency of the Method CN 
becomes quite obvious since the concentration of the tumor cells reaches 5.2 cells/mm but in Method IBE it is 2.25 cells/mm. 
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Table 1 

The typical behavior of the spatial relative error for At = 0.01 14.4 min. and At = 0.02 rs 28.8 min., 

Ax = 0.5 mm, c mM = 62.5 cells/mm, c™ 3 * = 39.89 cells/mm and p = 0.012 day -1 . 


Crank-Nicolson 


Implicit backward Euler 

Time (in days) 

/ 2 -norm 

Leo-norm 

/ 2 -norm 

Lqo- norm 


0.4334E—04 

0.5781E—04 

0.3760E—02 

0.5102E—02 

1 

0.8725E—06 

0.1819E—05 

0.3558E—03 

0.7649E—03 

5 

0.1621E—06 

0.4198E—06 

0.1330E—03 

0.3645E—03 

10 

0.4565E—08 

0.1132E—07 

0.1722E—04 

0.6324E—04 

50 

0.2694E—08 

0.4435E—08 

0.9755E—05 

0.3412E—04 

100 

0.4333E—08 

0.5946E—08 

0.1250E—04 

0.3121E—04 

150 

0.6307E—08 

0.753E—08 

0.1775E—04 

0.3488E—04 

200 

0.7873E—08 

0.9123E—08 

0.5548E—04 

0.5924E—04 

400 

0.3762E—08 

0.5284E—08 

0.6446E—04 

0.6884E—04 

500 

0.5713E—09 

0.8285E—09 

0.7235E—04 

0.8146E—04 

700 

0.1795E—09 

0.3235E—09 

0.5387E—04 

0.7320E—04 

800 

0.4273E—10 

0.1602E—09 

0.2557E—04 

0.3630E—04 

1000 

0.2569E—10 

0.1601E—09 

0.1605E—04 

0.2776E—04 

1050 

0.1503E—10 

0.1601E—09 

0.1047E—04 

0.206 IE-04 

1100 

0.1752E—10 

0.1601E—09 

0.5041E—05 

0.1053E—04 

1200 

0.0000E+00 

O.OOOOE+OO 

0.2575E—05 

0.5792E—05 

1280 


It took 114 min. with Method CN in 1280 days and with Method IBE in 2300 days to reach the maximum capacity c max = 
62.4989979 cells/mm, see Fig. 2(c) and (d). In Fig. 2(d), it is plotted only up to 1280 days. 

4.2. Stability condition 


It can be easily seen that Eq. (2.1) has two stationary states, namely c(t, x) = 0 and c(t, x) = c max . As usual it is assumed 
that the problem is posed on a bounded domain Bek with zero flux boundary conditions. Let us evaluate (in an ‘energy’ 
sense) 

d 
dt 


L 


| c(t, x)| 2 dx 


•■I 

■L 


3c(t, x) 

c(t, x)— 7 -dx 


c(t,x) 


dt 

3 ( dc(t, x)\ ( c(t, x) 

- D(x) + p C (t, x) 1 - | dx 

dx V dx ) \ Cmax 


Using integration by parts as usual and zero-flux conditions for the first integral on the right yields 


b 


— I |c(t,x)| 2 dx = -2 / D(x) 


■f [ 


dC(t, X) 

dx 


= 2 P [ 
Jb 


( c ( t ,*)) 2 


dx 
c(t,x) 


1 — 1 ) dx. 

Cmax 


Since the diffusion part gives a term with nonnegative sign, we have 


dt 


I, 


— / |c(t, x)| dx < 2 p / (c(t,x)r 1 - 


b 

i 


c(t,x) 


<2p (c(t, x)) 2 dx — 2p 


( c(t,x)Y 


fm 

Jb c, 


dx 


i 


< 2p / (c(t, x)) 2 dx 


(4.8) 


(4.9) 


(4.10) 


and therefore if 0 < c(0, x) < c max and p > 0, then f B |c(t, x)| 2 dx increases to the carrying capacity c max exponentially fast 
by the lemma on Gronwall’s inequality, see Figs. 4 and 7(a). 


4.3. Analysis 

We assume that the diffusion coefficient D depends on the tissue environment, thus we consider D as a function of x 
having finite number of discontinuities and this results in the lack of smoothness in our model solution. This is an obstacle 
for standard numerical procedures as we have explained in Section 3.1. 

The first two columns of Table 1 clearly show that in Method CN, l 2 and -norms of the spatial relative errors decrease 
much faster compared to the ones obtained by the Method IBE. Unless the stability condition is not violated, the spatial 
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Table 2 

The turning point for different carrying capacities and different initial concentrations, At = 0.01 
14.4 min. and p = 0.012 day -1 . 


Ax (mm) 

At (day) 

c mJX (cells/mm) 

Peak value of c(0, x) 

Time at turning point (day) 

0.5 

0.01 

62.5 

39.89 

42.7 

0.5 

0.01 

5000 

500 

41.87 

0.5 

0.01 

5000 

3000 

42.62 


relative errors for the Method IBE have the same decreasing pattern. Table 1 shows us that Method CN is not only more 
efficient in terms of errors but also much faster in terms of iteration steps. (This is because of its stability and accuracy.) 

Table 1 depicts the typical behavior of the spatial relative error in the case of Ax = 0.5 mm, p = 0.012 day -1 for 
At = 0.01 ~ 14.4 min. and At = 0.02 ~ 28.8 min. Similar results are obtained by using Ax = 0.25 mm. For fixed t 
values, the errors are calculated in different vector norms in order to form a comparison basis: 


and £/, 


Table 2 shows the turning points for each cases of the carrying capacities and the maximum initial tumor concentration 
of glioma cells. Taking At = 0.02 ~ 28.8 min also yields similar results. 

Fig. 4 shows that the tumor concentration of glioma cells, namely, c(t,x), reaches the maximum carrying capacity 
c max = 62.5 cells/mm which is the size of the stable, steady-state population when t = 1280 days by using the Method CN. 

Fig. 5(a) shows the maximum values of the tumor concentration of glioma cells, c^ ax (t) at each time step versus time for 
At = 0.01 ss 14.4 min., c max = 5000 cells/mm and c™ x = 3000 cells/mm and Fig. 5(b) is the zoomed graph of Fig. 5(a) 
focusing the turning point. As t increases, the tumor concentration of glioma cells reaches the maximum capacity that the 
tissue can bear. The case for c max = 62.5 cells/mm and c™ x = 39.89 cells/mm gives also a similar result graphically 
and numerically (an equivalent behavior observed). To verify the numerical methods, simulations have been performed for 
varying time-step, see Table 1 for Ax = 0.5 mm. 

To understand the effect of each term in our macroscopic reaction-diffusion evolution equation for cancer invasion, we 
plot the first term (net dispersal of glioma cells) and the second term (net proliferation of glioma cells) on the right-hand 
side ofEq. (2.1) versus time in days in Fig. 6 for each zone, respectively. It shows that until the day of 105 the net dispersal of 
glioma cells term is more dominant than net proliferation of glioma cells term in Zones 1 and 3, see Fig. 6(a) and (c), in other 
words, the tumor heterogeneity is even more noticeable as a result of increased abnormal cell proliferation. Between 105 
and 1000 days, net proliferation term dominates the net dispersal term. From Fig. 6, if we look at the data for the total net 
proliferation and net dispersal of glioma cells, we see that the rate of change of the tumor cell population density increases by 
t = 425 days and then it starts decreasing. This helps us to explain the convexity (concave up) and then concavity (concave 
down) occurring in the plot of the average number of tumor cells versus time, see Fig. 7. 


4.3.1. The average number of tumor cells, the mean radial distance, and the speed of the tumor cells 

To gain more insight of the tumor evolution, we also calculated the average number of tumor cells, mean radial distance, 
and the speed of the tumor cells for each zone and for the whole interval. Now, we define N(t), L(t), and S(t) which have 
the units of number of cells, the mean radial distance in mm, and the speed in mm/day, respectively, as following: the total 
number of tumor cells 

/»50 

N(t) := / c(t,x)dx, 

Jo 

the total mean radial distance of tumor cells 


L{t) 


/ 0 5 °xc(t, x) dx 

NO) 


[/ 50 x AffAi dx _ Ut ) J50 3c£x> d l 

and the total speed of tumor cells S(t) := = -L--1. 

We can see from Fig. 7(a) that Zone 2 (white region) contributes more to the number of cells since we assume that the 
tumor has grown to about 19.94 cells as a local mass before it begins to diffuse. Looking at the data plotted in Fig. 7(a), it 
tells us that in Zone 2, it took nearly 63 days to double the initial tumor cells. By t = 215 days, in Fig. 7(a), the tumor cells 
migrated further into the Zones 1 and 3 in the amount of the initial tumor cells in Zone 2. Additionally, by t = 218 days the 
average number of tumor cells in Zone 2 grew 10 times of the initial tumor cells and 20 times by t = 289 days. However, 
in Zone 1, the average number of tumor cells reached 398 cells by t = 622 days and in Zone 3, this was t = 574 days. In 
t = 631 days, the average number of cells grew 100 times of the initial tumor cells in Zone 2 (white region), 402 cells, 431 
cells for Zone 1 and Zone 3, respectively. 

In Fig. 7(b), by t = 27 days, the mean radial distance moves 10% to left in Zone 1 and to right 1.69% in Zone 3, just a little 
bit to left in Zone 2. As time evolves, by t = 103 days, in Zone 1 the mean radial distance moves 30% to left but 4.89% to right 
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(a) The peaks of concentration versus time. 



(b) The zoomed graph. 


Fig. 5. The maximum values of the tumor concentration cj^ft) of glioma cells at each time step versus time: p = 0.012 day At = 0.01 ss 14.4 min., 
c max = 5000 cells/mm, c™ ax = 3000 cells/mm and Ax = 0.5 mm using the Crank-Nicolson method. Fig. 5(b) shows the zoomed graph of Fig. 5(a) around 
turning point. 


Table 3 

The Crank-Nicolson Method: the infection days for each zone for different center position x 0 values of 
Gaussian initial profile (see Eq. (2.3)) for At = 0.01 14.4 min., Ax = 0.5 mm, c max = 62.5 cells/mm, 

c 0 max = 39.89 cells/mm and p = 0.012 day -1 . 


Gaussian profile position 
*o (mm) 

Infection day 



Zone 1 

Zone 2 

Zone 3 

(Zone 1)5 

- 

0.13 

37 

(Zone 2) 15 

1.29 

- 

17.41 

(Zone 2) 25 

7.69 

- 

6.81 

(Zone 2) 35 

18.83 

- 

0.98 

(Zone 3) 45 

39.16 

0.22 

- 


in Zone 3, see also Table 3. Moreover, in Fig. 7(b), the total mean radial distance graph shows that the tumor cells gather at 
the center of the Zone 2 (white matter), having the radius of 25 mm which means that tumor has occupied the whole tissue 
finally. 

Even though the speed of tumor cells changes sign for each zone, the total speed of tumor cells does not change, see 
Fig. 7(c). After t = 800 days, the total number of tumor cells become constant («3090) cells as shown in Fig. 7(a), thus the 
total speed of tumor cells is zero, see Fig. 7(c). 

In all the computational runs, we considered c max = 62.5 cells/mm and took the initial condition as a Gaussian initial 

, iflzlol 2 

tumor profile c(0, x) = -r^-e 2 v 8 / where x 0 = 25 mm (the middle of the considered interval) and e = 0.01. Table 3 
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1.4 


7 


1.2 

1 

0.8 

0.6 

0.4 

0.2 



-Net Dispersal Term 

Net Proliferation Term 


6 

5 

4 

3 

2 

1 

0 


,' \ Net Dispersal Term 

/ \ Net Proliferation Term 



0 L ——- 1 - 1 -- 1 — — ~ - - 1 - -1 - 1 - 1 - 1 - 1 - 1 - 1 - 

0 200 400 600 800 1000 1200 1400 0 200 400 600 800 1000 1200 1400 


Time (days) 

(a) Zone 1 (gray region): 0 < x < 7.5 mm. 


Time (days) 

(b) Zone 2 (white region): 7.5 < x < 42.5 mm. 


1 


0.5 


0 

0 200 400 600 800 1000 1200 1400 

Time (days) 

(c) Zone 3 (gray region): 42.5 < x < 50 mm. 

Fig. 6. The net dispersal term (DT) and the net proliferation term (PT) of glioma cells versus time in days of Eq. (2.1) for p = 0.012 day -1 , At = 0.01 « 
14.4 min., c max = 62.5 cells/mm, c™ ax = 39.89 cells/mm and Ax = 0.5 mm. 

shows the infected days for each zone for different center position numbers x 0 ranging from 0 to 50. If the center position 
number x 0 = 5, i.e. the center of Gaussian initial profile is located in Zone 1, then the average number of cells for Zone 1 is 
larger than the one in Zone 3, and the net dispersal term in Zone 1 is more than the one in Zone 3, but as time gets larger 
they reach an equilibrium value. These observations can be reversed if the center of Gaussian initial profile is located in Zone 
3 rather than in Zone 1. 

5. Conclusion 

In this paper, a mathematical model using one dimensional parabolic equation regarding the growth of human tumors has 
been presented. What makes this study different from the existing literature is the fact that we check the stability condition 
and the accuracies of the numerical methods. Since the diffusion coefficient D has a finite number of discontinuities, our 
model solution has also a lack of smoothness. However, at the steady-state, these discontinuities vanish. Two finite difference 
methods are employed, namely, the implicit backward Euler method and the Crank-Nicolson method. The Crank-Nicolson 
implicit method is shown to be the better choice to solve the equation. 

Our model shows that tumor cell population reaches a maximum cell number that the tissue can carry, however, in reality, 
as the cells reach this “steady state", new mutations will continue to occur within the nuclear DNA. These modifications will 
lead to new phenotypic features of the cancer cell, such as having higher proliferation rates or being more invasive. 


Net Dispersal Term 
Net Proliferation Term 


/ \ 

/ \ 
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15 


-Total mean radial distance of tumor cells 

Zone 1 (Grey region) 

--Zone 2 (White region) 

- - Zone 3 (Grey region) 


10 - 

5 - . 

0 - 1 - 1 - 1 - 1 - 1 - 1 - 
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Time (days) 


Q 0.02 


0.015 

0.01 

0.005 
S(t) 

0 

-0.005 

- 0.01 

-0.015 

0 200 400 600 800 1000 1200 1400 

Time (days) 

Fig. 7. For each zone: (a) the average number, N(t), of tumor cells, (b) the mean radial distance, L(t), of tumor cells, (c) the speed, S(t), of tumor cells. All 
these calculations are done for p = 0.012 day -1 , At = 0.01 -- 14.4 min., c max = 62.5 cells/mm, c™ 3 * = 39.89 cells/mm and Ax = 0.5 mm. 

We have also numerically studied the behavior of the evolution of tumor concentration of the glioma in terms of the 
number of tumor cells, the average mean radial distance and the speed of tumor cells for different center position values of 
Gaussian initial profile for each zone. 

For future work, we plan to solve the problem in three dimensions by including the radiotherapy and chemotherapy 
factors. This will help developing ways to integrate data on a patients cancer into personalized models of tumor growth. 
Such models can better predict how a tumor will respond to different treatments and drugs than any one piece of data 
alone. 
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